function see_igvs(v_vec,s_vec,flag_spec,lb,ub)
% See the IG distribution with parameters (v,s) 
% where flag_spec='L' if using Lubrano's notation and 
%       flag_spec='J' if using Judge's   notation 
% Note that s(L)=v*s(J)^2
if nargin < 3 
    flag_spec='L'; 
    disp('Default is Lubrano et al. specification') 
end 
if nargin < 4 
    lb=0.001; 
    ub=3; 
end 
if strcmp( upper(flag_spec),'J' )~=1 && strcmp( upper(flag_spec),'L' )~=1 
    error('Must choose either L(ubrano) or J(udge) specification') 
end 
nrow=length(v_vec); 
if length(s_vec)~=nrow;error('Size of V_VEC and S_VEC must match');end 

npoints=10000; 
grid=linspace(lb,ub,npoints); 
momat=zeros(3,nrow); 
pdfmat=zeros(npoints,nrow); 

figure; 

lvec={'-','--','.-'}; 
wvec=[1.5;2;1.5]; 

for ii=1:nrow; 
    v=v_vec(ii); 
    s=s_vec(ii); 
    dispaj('Begin V=',v,' , S=',s); 
    
    cons=2/(gamma(v/2)); 
    if strcmp( upper(flag_spec),'J' )
        cons=cons*( (v*s*s/2)^(0.5*v) );
        expcons=-0.5*(v*s*s);
        m1=(s*sqrt(v/2)*gamma((v-1)/2))/( gamma(v/2) );
        m2=sqrt((v*s*s*(m1*m1))/(v-2));
        maxv=s*((v/(v+1))^(1/2));
        for jj=1:npoints; 
            pdfmat(jj,ii)=( (grid(jj))^( -(v+1) ) )*exp(expcons/(grid(jj)^2)); 
        end 
        pdfmat(:,ii)=pdfmat(:,ii)*cons; 
    end 
    momat(1,ii)=m1; 
    momat(2,ii)=m2; 
    momat(3,ii)=maxv; 
    
    dispaj('Mean=',m1,'    ','Std=',m2,'   Mode=',maxv); 
    
    plot(grid(:),pdfmat(:,ii),lvec{ii},'LineWidth',wvec(ii)); 
    if ii==1; hold on; end 
       
end 
hold off; 
            